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1. Introduction 

In the past decade the progress of simulation algorithms and the availability of more affluent 
computational resources have enabled us to investigate two- and three-flavor dynamical quark ef- 
fects with lighter quark masses. The CP-PACS and JLQCD joint project[|l]] has performed a 2+1 
flavor full QCD simulation employing the 0(a) -improved Wilson quark action with nonperturba- 
tive csw and the Iwasaki gauge action^] on a (2fm) 3 box at three lattice spacings. Unfortunately, 
the lightest up and down quark mass reached is about 64MeV corresponding to m n /m p w 0.6. 

A goal of the next decade should be the realistic full QCD simulation with the physical up, 
down and strange quark masses. The PACS-CS project, which is based on the PACS-CS (Parallel 
Array Computer System for Computational Sciences) computer installed at University of Tsukuba 
on July 1 of 2006 P], aims at this goal succeeding the Wilson-clover 2+1 program of the previous 
CP-PACS/JLQCD project. 

In this report we present a preparatory study for the PACS-CS project using (1.6fm) 3 lattice. 
To simulate the dynamical up and down quarks, we employ the domain-decomposed HMC algo- 
rithm proposed by Liischer[|]]. The effectiveness of this algorithm for small quark mass region is 
already shown in the two-flavor case[Q, ||]. The strange quark is included by the exact Polynomial 
HMC (PHMC) algorithm Our primary purpose is to investigate how light up and down quark 
masses we can go down to with LDDHMC. The efficiency of LDDHMC is compared to that of 
the conventional HMC. Based on this study we finally make a parameter choice for the production 
runs on PACS-CS, and briefly discuss the physics plan. 

2. Simulation details 

We employ the 0(a) -improved Wilson quark action with nonperturbative csw[0] and the 
Iwasaki gauge action at j8 = 1.9 on a 16 3 x 32 lattice. The lattice spacing is about O.lfm from 
the previous CP-PACS/JLQCD results on a 20 3 x 40 lattice[JT|]. In Table @ we summarize the de- 
tails of the simulation parameters. We employ four degenerate up and down quark masses based on 
the previous CP-PACS/JLQCD results, while the strange quark mass is fixed at fc s = 0.1364 close 
to the physical point K s = 0.136412(50)[j|]. The LDDHMC algorithm is implemented with a 8 4 
block size, while the exact PHMC algorithm is not domain-decomposed. We choose the trajectory 
length of t = 0.5 /y/2. With the integers No, N\, N2 the step sizes are given by Sx = xj (N0N1N2) 
for the gauge part, 8x = x/(N\N2) for the strange quark and the UV part of the up and down 
quarks, and 8x = X /N2 for the IR part of the up and down quarks. The Wilson quark matrix in- 
version is carried out by the even/odd preconditioned BiCGStab solver with the stopping condition 
\Dx — b\/\b\ < 10~ 9 for the force calculation and 10~ 14 for the Hamiltonian. 

All the simulations presented in this report are carried out on Hitachi SR11000/J1 at Informa- 
tion Technology Center of the University of Tokyo. 

3. Numerical results 
3.1 History of dH 

In Fig.|I| we show the history of dH for each fc^. We observe a frequent occurrence of spikes 
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Table 1: Simulation parameters and basic physical quantities. #mult is the number of multiplications of the 
Wilson-Dirac quark matrix on the full lattice including even and odd sites. 
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for fc u( j = 0.13759 and 0.13770, whose probabilities for \dH\ > 2 are about 2%. Although the 
acceptance ratios for k u< i = 0.13759 and 0.13770 are still sufficiently high, {t~ dH ) which measures 
the area preserving property deviates from unity by 2a. It is clear that the spikes are a potential 
source to violate the area preserving property. In order to avoid any problem associated with the 
spikes, we shall incorporate the replay trick[^, ||] in the production run. 

3.2 Magnitude of force terms 

Figure |2| shows the magnitude of the four force terms as a function of the axial vector Ward 
identity quark mass discussed below; F g denotes the gauge part, Fuv,ir for the UV and IR 
parts of the up and down quarks, and F s for the strange quark. Note that the definition of the 
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Figure 1: Time histories of dH as a function of the trajectory number after thermalization for fc u d 
0.13700(topleft), 0.13741(bottom left), 0.1 3759(top right) and 0.1 3770(bottom right). 



magnitude of force term is a factor two larger than that in Ref.[Q]. We observe that \\F g \\, H^uvll 
and \\F S \\ are almost independent of the quark mass, while ||Fir|| gradually increases as the up and 
down quark masses decreases. At K a & = 0.13770 we find 



\F g \\ : H-Puvl 



|F s ||:||Fir||«24:6:3:1. (3.1) 
This result suggests that our choice of 8x for the strange quark may be unnecessarily fine. 

3.3 Comparison of LDDHMC and HMC 

For a direct comparison between LDDHMC and HMC we have repeated the 2+1 flavor sim- 
ulation at k u & = 0.13700 and fc s = 0.13640 with the conventional HMC algorithm employing the 
same parameters except T = 1 for HMC. Analyzing 3000 trajectories after thermalization we obtain 
Tint [P] = 6.6(1.8) and #mult/traj= 148625(626) for HMC. Since the trajectory lengths are differ- 
ently chosen for LDDHMC and HMC, we compare their efficiency by v mi [P] -#mult/10 3 which 
should be independent of the trajectory length. This quantity is found to be 215(55) for LD- 
DHMC and 981(268) for HMC. Hence LDDHMC is almost five times more efficient than HMC 
at mps/my ~ 0.6. We should keep in mind that this comparison involves the strange quark part 
simulated by the exact PHMC algorithm. 

3.4 Quark masses 

The up and down quark masses are measured by using the axial vector Ward identity (AWI): 

AWI _ Z A <0|V4A; mpr |PS> 



ud 



Zp 2{0\P\PS) 



(3.2) 
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Figure 2: Magnitude of each force term as a function of the AWI quark mass m^ 1 . The errors are within 
symbols. 



where the renormalization factors Za,p and the improvement coefficients are determined pertur- 
batively up to one-loop level. The results are given in Table [T] in physical units together with 
the corresponding pseudoscalar masses. It is encouraging that LDDHMC allows a simulation at 
ra™ = 15MeV which is roughly a quarter of the lightest quark mass employed in the previous 
CP-PACS/JLQCD project. In Table @ we also give the results for nip S /m^ 1 , whose small quark 
mass dependence indicates that the finite volume effects are not sizable. 

3.5 Spectral gap distribution 

Recently Del Debbio et al. have discussed the relevance of the spectral gap distribution of 
the Wilson-Dirac operator for the stability of two-flavor lattice QCD simulations^. This study is 
applied to our three flavor case. We are interested in the effects of the strange quark contributions. 

Following Ref.[]9|] we use the hermitian operator Q m = JsD m rather than the Wilson-Dirac 
operator D m = (l/2){y fl (V* + V^) — aV* V^} + niQ. The spectral gap is defined as 

H = min{|A| | A is an eigenvalue of Q m }. (3.3) 

In Fig.|| we show histograms of the spectral gap /I for our four values of fc u d, which are obtained 



by the implicitly restarted Lanczos algorithm[10]. We observe roughly symmetric distributions for 



all the quark masses. Their median and average given in Table [j] are consistent with each other 
within error bars. In Fig.|] we plot the quark mass dependencies for /Vedian an d /^average- Although 
both quantities are roughly proportional to the AWI quark mass, we find a clear tendency that they 
deviate upwards from the linearity in terms of m^ 1 toward the chiral limit. 

For the width of the distribution a we employ the same definition as in Ref.[]9|]: a is defined as 
(v — u)/2, where [u,v] is the smallest range of pL which contains more than 68.3% of the data. This 
is to avoid potentially large statistical uncertainties which might occur when data are not sufficiently 
sampled. The results of a are given in Table [jj We find that the width of the distribution diminishes 
as the up and down quark mass decreases. An intriguing quantity is the combination Gy/V /a: In 
the two-flavor case Del Debbio et al. found that its value is roughly unity independent of the quark 
mass, volume, and lattice spacingTOn. In our three-flavor case, on the contrary, we find GW/a < 1 
and its value decreases as the up and down quark mass decreases. This fact may suggest that the 
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Figure 3: Histograms of fx for K ud = 0.13700(left top), 0.13741(left bottom), 0.13759(right top) and 
0.13770(right bottom). Vertical dashed lines denote jU me dian- 



contributions of the strange quark stabilize the simulation at the lighter quark masses by shrinking 
the gap distribution of the Wilson-Dirac operator. 

4. Plan for production run on PACS-CS 

We employ the LDDHMC algorithm with the replay trick for the up and down quarks. The 



Wilson-Dirac quark matrix inversion on the full lattice is accelerated by the SAP+GCR solver] 11], 



which is three times faster than the BiCGStab algorithm for our parameter range[|ll|]. For the 



strange quark we use the UV-filtered PHMC algorithm, which is two to three times more efficient 



than the PHMC algorithm[[12|] 



We choose j8 = 1.83, 1.90, 2.05 which are the same as in the previous CP-PACS/JLQCD 
project[|]. The physical lattice volume is fixed at (3.0fm) 3 : 24 3 x 48 for j8 = 1.83, 32 3 x 64 for 
j8 = 1.90 and 40 3 x 80 for j8 = 2.05. Assuming that the relation a °c a/s/V found in Ref.[g] for 
the two-flavor simulation is applicable to the three-flavor case, we expect that a stable simulation at 
m™ = 4MeV is possible at j3 = 1 .90 and even smaller up and down quark masses are accessible at 
/3 = 2.05. Two strange quark masses are employed to interpolate the results on the physical point. 
We plan to accumulate 100 independent configurations with 10 4 trajectories. 

The first step of the project is a simulation with the up and down quark mass down to m^* 1 = 
7MeV, which is about twice the physical mass and one order smaller than the lightest quark mass 
employed in the previous CP-PACS/JLQCD project. Once this is achieved, we can almost remove 
the systematic errors associated with the chiral extrapolation which are currently one of the most 
significant errors in lattice QCD simulations. The hadron spectrum together with the quark masses 
should be determined precisely. We are now preparing a nonperturbative renormalization factor for 
the quark masses employing the Schrodinger functional method. 
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Figure 4: Median and average of the gap distribution as a function of AWI quark mass fflCj . The dashed 
line drawn from the origin to ^median at the largest quark mass is to guide the eye. 



There are two physical subjects which are especially interesting for the light quark masses. 
One is physics associated with topology. The CP-PACS/JLQCD Collaborations presented an en- 
couraging result for the r\' meson mass[[l3|]. An intriguing point is that the signals become cleaner 



as the up and down quark mass decreases. The other is the 1 = K-K scattering. In Ref.[|14]] we 
discussed that the light dynamical quarks are essential ingredients to study it. 

This work is supported in part by Grants-in-Aid for Scientific Research from the Ministry 
of Education, Culture, Sports, Science and Technology (Nos. 13135204, 13135216, 15540251, 
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